#include <stdlib.h>
#include <stdio.h>
#include <omp.h>

void mv_orig(int m, int n, double M[m][n], double V[n], double W[m])
{
  int i, j;

  for (i = 0; i < m; i++) {
    W[i] = 0.0;
    for (j = 0; j < n; j++) {
      W[i] += M[i][j] * V[j];
    }
  }
}

double reduce(int n, double M_i[n], double V[n], int actual);
void mv(int m, int n, double M[m][n], double V[n], double W[m])
{
  int i, j;

  for (i = 0; i < m; i++) {
    W[i] = reduce(n, M[i], V, 0);
  }
}

double reduce(int n, double* M_i, double* V, int actual) {
    // recursion tail, if indices coincide computation is trivial
    if (n == 1) {
        return M_i[0] * V[0];
    }

    // halve workload using a recursive call, with parallel for O(log n) time steps
    double sum[2];
    for (int i = 0; i < 2; i++) {
    sum[i] = reduce(
            n / 2 + (i * (n % 2)),
            &M_i[(n / 2) * i],
            &V[(n / 2) * i],
            actual + (n / 2) * i
    );
  }

  return sum[0] + sum[1];
}

void funky_sort(int n, int a[n]) {
    int b[n];

    #pragma omp parallel for shared (a, b)
    for (int i = 0; i < n; i++) {
        int j;
        int count = 0;

        for (j = 0; j < i; j++) {
            if (a[j] <= a[i]) count++;
        }
        j++;
        for (; j < n; j++) {
            if (a[j] < a[i]) count++;
        }
        b[count] = a[i];
    }

    #pragma omp parallel for shared (a, b)
    for (int i=0; i<n; i++)
        a[i] = b[i];
}

//int main(int argc, char *argv[]) {
int main(void) {
    int max_threads = omp_get_max_threads();

    printf("max number of threads is %d\n", max_threads);

    #pragma omp parallel for
    for (int i = 0; i < max_threads * 2; i++) {
        printf("hello from thread %d\n", omp_get_thread_num());
    }

    #pragma omp parallel
    {
        printf("hello %d\n", omp_get_thread_num());
        printf("again %d\n", omp_get_thread_num());
    }

    double M[9][9] = {
        { 1, 2, 3, 4, 5, 6, 7, 8, 9 },
        { 1, 2, 3, 4, 5, 6, 7, 8, 9 },
        { 1, 2, 3, 4, 5, 6, 7, 8, 9 },
        { 1, 2, 3, 4, 5, 6, 7, 8, 9 },
        { 1, 2, 3, 4, 5, 6, 7, 8, 9 },
        { 1, 2, 3, 4, 5, 6, 7, 8, 9 },
        { 1, 2, 3, 4, 5, 6, 7, 8, 9 },
        { 1, 2, 3, 4, 5, 6, 7, 8, 9 },
        { 1, 2, 3, 4, 5, 6, 7, 8, 9 },
    };

    double V[9] = { 1, 2, 3, 4, 5, 6, 7, 8, 9 };
    double W[9];
    double W2[9];
    mv(9, 9, M, V, W);
    mv_orig(9, 9, M, V, W2);
    for (int i = 0; i < 9; i++) {
        printf("%d %f %f\n", i, W[i], W2[i]);
    }

    int funky_sort_test[9] = { 123, 1328, 132, 19325, 538, 13589, 248, 24187, 248 };
    funky_sort(9, funky_sort_test);
    for (int i = 0; i < 9; i++)
        printf("%d %d\n", i, funky_sort_test[i]);

    return 0;
}
